#' Get MODIS data and process to match geogrid
#' 
#' \code{GetMODIS} downloads, mosaics, and resamples MODIS data to match input 
#' geogrid.
#' 
#' \code{GetMODIS} reads a geogrid file and parameters on MODIS product and data
#' range and downloads MODIS tiles where they do not exist (using the R 
#' \href{http://r-forge.r-project.org/projects/modis/}{MODIS} package), mosaics 
#' where necessary, clips, reprojects, and resamples (using nearest neighbor) to
#' match the geogrid. Results in a set of TIF files per the 
#' \href{http://r-forge.r-project.org/projects/modis/}{MODIS} package 
#' specifications.
#' 
#' Please see documentation on the R 
#' \href{http://r-forge.r-project.org/projects/modis/}{MODIS} package for 
#' details on required installs and workspace setup. This tool builds off of the
#' \href{http://r-forge.r-project.org/projects/modis/}{MODIS} 
#' \code{\link[MODIS]{runGdal}} tool, so follows 
#' \href{http://r-forge.r-project.org/projects/modis/}{MODIS} file directory 
#' structure (see \code{\link[MODIS]{MODISoptions}}). This tool requires a local
#' \href{http://www.gdal.org/}{GDAL} installation in addition to the required R 
#' packages.
#' 
#' NOTE: This tool currently only works for geogrid files in Lambert Conformal 
#' Conic projection OR a GDAL-friendly georeferenced TIF file.
#' 
#' @param extent The pathname to the geogrid file (i.e., geo_em.d01.nc)
#'  OR the path to a georeferenced TIF file of the domain OR the path to 
#'  a shapefile of the desired extent OR a standard country name that can
#'  be identified in the mapdata package.
#' @param prodName The MODIS product name to download/process. Run the 
#'   \href{http://r-forge.r-project.org/projects/modis/}{MODIS} package 
#'   getProducts() for a complete list of supported products.
#' @param outDir Directory name to store processed TIF files. This is the 
#'   equivalent to the 
#'   \href{http://r-forge.r-project.org/projects/modis/}{MODIS} package's "job" 
#'   name. This is a directory name only and NOT a full path. The directory will
#'   be created in the preset 
#'   \href{http://r-forge.r-project.org/projects/modis/}{MODIS} package 
#'   outDirPath.
#' @param begin Date string for the start date to download/process MODIS tiles. 
#'   The date string should follow 
#'   \href{http://r-forge.r-project.org/projects/modis/}{MODIS} package 
#'   convention (e.g., "2011.06.01").
#' @param end Date string for the end date to download/process MODIS tiles. The 
#'   date string should follow 
#'   \href{http://r-forge.r-project.org/projects/modis/}{MODIS} package 
#'   convention (e.g., "2011.06.01").
#' @param collection From MODIS::runGdal help: "Default is to download most
#'   recent collection version. See: ?getCollection"
#' @param buffer From MODIS::runGdal help: "Numeric. Buffer [in units of the
#'   'outProj'] around the specified extent. See: ?getTile" (DEFAULT=0.04, which
#'   seems to cover MODIS tile corners sufficiently)
#' @param SDSstring From MODIS::runGdal help: "Default is extract all SDS
#'   (layers). See: ?getSds."
#' @param exclList List pairing SDS name (e.g., "Lai_1km", "FparLai_QC") with 
#'   values to exclude from interpolation (i.e., convert to "nodata"). For 
#'   example, for the Lai_1km product, valid range is 0-100, values 249-255 are 
#'   fill values, but only one value (255) is reported as _FillValue in the HDF 
#'   metadata. The exclude list should specify value ranges to exclude, and 
#'   should be in the form: \deqn{list("SDS name"="exclude range", ...)} where 
#'   "exclude range" can be a single number (e.g., 254), "gt <somenumber>" to 
#'   exclude all values over some threshold (e.g., "gt 100"), or "lt 
#'   <somenumber>" to exclude all values less than some threshold (e.g., "lt 
#'   0").
#' @param resampList List pairing SDS name (e.g., "Lai_1km", "FparLai_QC") with 
#'   reasmpling method to use. GDAL possible options (depends on version): 
#'   'near', 'bilinear', 'cubic', 'cubicspline', 'lanczos', 'mode', 'average' 
#'   (DEFAULT='near', unless otherwise specified in 
#'   \code{MODISoptions(resamplingType="<somemethod>")}). The resampling list 
#'   should be in the form: \deqn{list("SDS name"="method", ...)} where "method"
#'   is one of the GDAL methods listed above.
#' @param quiet From MODIS::runGdal help: Logical, Default FALSE. Some progress 
#'   informations.
#' @param scriptPath OPTIONAL path to where gdal_calc.py script lives in case
#'   it is not in the same default path as gdal executables
#' @return Empty
#'   
#' @examples
#' ## First, specify the target directories for the MODIS package.
#' \dontrun{
#' MODISoptions(localArcPath="/d1/WRF_Hydro/RS/MODIS_ARC", 
#'              outDirPath="/d1/WRF_Hydro/RS/MODIS_ARC/PROCESSED",
#'              stubbornness="low")
#' 
#' ## Then, use GetMODIS to download the MODIS MOD15A2 (FPAR/LAI) product for 
#' ## all tiles that overlap our Fourmile Creek domain (as specified by the 
#' ## geo_em.d01.nc file), mosaic if multiple tiles, clip to the domain extent,
#' ## and resample to the domain grid cells. The raw HDF files will be stored in
#' ## /d1/WRF_Hydro/RS/MODIS_ARC/ and the final processed TIF files will be 
#' ## stored in /d1/WRF_Hydro/RS/MODIS_ARC/PROCESSED/Fourmile_LAI/.
#' 
#' GetMODIS(extent="/d1/WRF_Hydro/Fourmile_fire/DOMAIN/geo_em.d01.nc", 
#'          prodName="MOD15A2", outDir="Fourmile_LAI", 
#'          begin="2011.01.01", end="2011.01.31",
#'          exclList=list("Fpar_1km"="gt 100", 
#'                        "Lai_1km"="gt 100", 
#'                        "FparLai_QC"="255", 
#'                        "FparExtra_QC"="255", 
#'                        "FparStdDev_1km"="gt 100", 
#'                        "LaiStdDev_1km"="gt 100"), 
#'          resampList=list("Fpar_1km"="bilinear", 
#'                        "Lai_1km"="bilinear", 
#'                        "FparLai_QC"="mode", 
#'                        "FparExtra_QC"="mode", 
#'                        "FparStdDev_1km"="bilinear", 
#'                        "LaiStdDev_1km"="bilinear"))
#' }
#' @keywords IO
#' @concept MODIS dataGet
#' @family MODIS
#' @export
GetMODIS <- function(extent, prodName, outDir, begin=NULL, end=NULL, 
                     collection=NULL, buffer=0.04,
                     SDSstring=NULL, 
                     exclList=NULL, resampList=NULL, quiet=FALSE, scriptPath=NULL) {
    # Check packages
    if (!(require("rgdal") & require("raster") & require("ncdf4") & require("MODIS"))) {
        stop("Required packages not found. Must have R packages: rgdal (requires GDAL system install), raster, ncdf4, and MODIS")
        }
    # Get paths
    locPath <- paste0(options("MODIS_localArcPath"))
    if (!file.exists(locPath)) {
      dir.create(locPath, showWarnings = FALSE)
    }
    # Get geogrid and projection info
    if (inherits(try(suppressWarnings(ncdf4::nc_open(extent))), "try-error")) {
        if (file.exists(extent)) {
            if ( grepl(".shp", extent, ignore.case = TRUE) ) {
                message("Assuming file is a shapefile.")
                hgt.r <- extent
            } else if ( grepl(".tif", extent, ignore.case = TRUE) ) {
                message("Assuming file is a georeferenced TIF.")
                hgt.r <- raster::raster(extent)
            } else {
                stop("Must provide a .shp or .tif file.")
            }
        } else {
            message("Assuming file is a country name. Requires 'mapdata' package.")
            hgt.r <- extent
        }
    } else {
      message("Assuming file is a geogrid in netcdf format.")
      geogrd.nc <- ncdf4::nc_open(extent)
      map_proj <- ncdf4::ncatt_get(geogrd.nc, varid=0, attname="MAP_PROJ")$value
      cen_lat <- ncdf4::ncatt_get(geogrd.nc, varid=0, attname="CEN_LAT")$value
      cen_lon <- ncdf4::ncatt_get(geogrd.nc, varid=0, attname="STAND_LON")$value
      truelat1 <- ncdf4::ncatt_get(geogrd.nc, varid=0, attname="TRUELAT1")$value
      truelat2 <- ncdf4::ncatt_get(geogrd.nc, varid=0, attname="TRUELAT2")$value
      if (map_proj==1) {
         geogrd.crs <- paste0("+proj=lcc +lat_1=", truelat1, " +lat_2=", truelat2, 
                              " +lat_0=", cen_lat, " +lon_0=", cen_lon, 
                              " +x_0=0 +y_0=0 +a=6370000 +b=6370000 +units=m +no_defs")
        } else {
        stop('Error: Projection type not supported (currently this tool only works for Lambert Conformal Conic projections).')
        }
    dx <- ncdf4::ncatt_get(geogrd.nc, varid=0, attname="DX")$value
    dy <- ncdf4::ncatt_get(geogrd.nc, varid=0, attname="DY")$value
    if ( dx != dy ) {
        stop(paste0('Error: Asymmetric grid cells not supported. DX=', dx, ', DY=', dy))
        }
    # Create a readable TIF from geogrid
    ExportGeogrid(extent, "HGT_M", paste0(locPath, "/geogrid_tmp.tif"))
	  hgt.r <- raster::raster(paste0(locPath, "/geogrid_tmp.tif"))
    system(paste0("rm ", paste0(locPath, "/geogrid_tmp.tif")))
    }
    # Run the download & processing
    mod.list <- MODIS::runGdal(product=prodName, collection=collection, 
                               begin=begin, end=end,
                               extent=hgt.r, buffer=buffer, 
                               SDSstring=SDSstring, job=outDir, 
                               exclList=exclList, resampList=resampList,
                               quiet=quiet, scriptPath=scriptPath)
}



#' Convert a set of MODIS images to a raster stack and, optionally, a NetCDF
#' file.
#' 
#' \code{ConvertRS2Stack} takes a set of pre-processed MODIS TIF files and
#' creates a raster stack and, optionally, a NetCDF file.
#' 
#' \code{ConvertRS2Stack} scans the specified directory and imports
#' pre-processed MODIS TIF files matching the specified expression and combines
#' them into an R raster stack object and, optionally, an output NetCDF file.
#' Files in the input directory should be already processed through the
#' \code{\link{GetMODIS}} tool or follow the same file naming convention used by
#' the MODIS \code{\link[MODIS]{runGdal}} tool. See the R
#' \href{http://r-forge.r-project.org/projects/modis/}{MODIS} package for more
#' on those specifications.
#' 
#' @param inPath The path to the directory that holds the already processed
#'   MODIS TIF files.
#' @param matchStr The regular expression for filename matching (e.g.,
#'   "*Lai_1km.tif").
#' @param begin Date string for the start date to include. The date string
#'   should follow \href{http://r-forge.r-project.org/projects/modis/}{MODIS}
#'   package convention (e.g., "2011.06.01"). (DEFAULT=NULL, all files are
#'   processed).
#' @param end Date string for the end date to include. The date string should
#'   follow \href{http://r-forge.r-project.org/projects/modis/}{MODIS} package
#'   convention (e.g., "2011.06.01"). (DEFAULT=NULL, all files are processed).
#' @param noData Value for specifying "no data" per the MODIS product. Should be
#'   combined with a qualifier. For example, for the MOD15A2 (LAI) product,
#'   valid value range is 0-100. So the noData value should be set to 100 and
#'   the noDataQual to "max". (DEFAULT=NULL, no additional "no data" screening
#'   is applied)
#' @param noDataQual Qualifier for specifying "no data" per the MODIS product.
#'   Should be combined with a nodata value. For example, for the MOD15A2 LAI
#'   product, valid value range is 0-100. So the noData value should be set to
#'   100 and the noDataQual to "max". Options are: "exact" (value == noData are
#'   converted to NA), "min" (value < noData are converted to NA), "max" (value
#'   > noData are converted to NA). (DEFAULT="exact", only exact matches are
#'   converted to NA)
#' @param valScale The scale factor to apply to the image values per the MODIS
#'   product. The adjustment is applied as VALUE * valScale + valAdd. For
#'   example, for the MOD15A2 LAI product, the scale factor is 0.1. (DEFAULT =
#'   1)
#' @param valAdd The addition value to apply to the image values per the MODIS
#'   product. The adjustment is applied as VALUE * valScale + valAdd. (DEFAULT =
#'   0)
#' @param outFile OPTIONAL name for an output NetCDF file. A NetCDF file will
#'   only be created if this file name is provided. Images will be exported
#'   after the "no data" and scale adjustments are made. If you want to do
#'   smoothing or other time series processing, do not export a NetCDF file here
#'   but do processing on the raster stack and then use ConvertStack2NC to
#'   export the processed brick to a NetCDF file.
#' @param varName Name for the NetCDF export variable. Only required if outFile
#'   is provided.
#' @param varUnit Units for the NetCDF export variable. Only required if outFile
#'   is provided.
#' @param varLong Long name for the NetCDF export variable. Only required if
#'   outFile is provided.
#' @param varNA Value to set for "NA" or "no data". Default is -1.e+36. Only
#'   required if outFile is provided.
#' @param pos1 From MODIS orgTime: Start position of date in the filename 
#'   (DEFAULT=10).
#' @param pos2 From MODIS orgTime: End position of date in the filename 
#'   (DEFAULT=16).
#' @param format From MODIS orgTime: How is the date formatted in the file.
#'   Default is "\%Y\%j" ('YYYYDDD'). Read ?as.Date for for more information.
#' @return A raster stack.
#' @examples
#' ## Import the already processed LAI TIF images into a raster stack. Use 
#' ## the full time series of images.
#' 
#' \dontrun{
#' lai.b <- ConvertRS2Stack("/d6/adugger/WRF_Hydro/RS/MODIS_ARC/PROCESSED/BCNED_LAI",
#'                          "*Lai_1km.tif", noData=100, noDataQual="max", 
#'                          valScale=0.1, valAdd=0)
#' 
#' ## Export a subset of the already processed LAI TIF images into an output netcdf file
#' 
#' lai.b <- ConvertRS2Stack("/d6/adugger/WRF_Hydro/RS/MODIS_ARC/PROCESSED/BCNED_LAI", 
#'                          "*Lai_1km.tif", begin="2011.06.01", end="2011.06.30", 
#'                          noData=100, noDataQual="max", valScale=0.1, valAdd=0, 
#'                          outFile="BCNED_LAI.nc", varName="LAI",
#'                          varUnit="(m^2)/(m^2)", varLong="Leaf area index")
#' }
#' @keywords IO
#' @concept MODIS dataMgmt
#' @family MODIS
#' @export
ConvertRS2Stack <- function(inPath, matchStr, begin=NULL, end=NULL,
                            noData=NULL, noDataQual="exact", valScale=1, valAdd=0,
                            outFile=NULL, varName=NULL, varUnit=NULL, varLong=NULL, 
                            varNA=(-1.e+36),
                            pos1=10, pos2=16, format="%Y%j") {
    # Get file list
    if (!is.null(begin) & !is.null(end)) {
        beginDt <- format(as.POSIXct(begin, format="%Y.%m.%d", tz="UTC"), 
                          "%Y%j", tz="UTC")
        endDt <- format(as.POSIXct(end, format="%Y.%m.%d", tz="UTC"), 
                        "%Y%j", tz="UTC")
        timeInfo <- MODIS::orgTime(MODIS::preStack(path=inPath, pattern=matchStr), 
                                   begin=beginDt, end=endDt, pillow=0,
                                   pos1=pos1, pos2=pos2, format=format)
    } else {
        timeInfo <- MODIS::orgTime(MODIS::preStack(path=inPath, pattern=matchStr), 
                                   pillow=0, pos1=pos1, pos2=pos2, format=format)
    }
    rsFilelist <- MODIS::preStack(path=inPath, pattern=matchStr, timeInfo=timeInfo)
    #rsFilelist <- list.files(path=inPath, pattern=glob2rx(matchStr), full.names=TRUE)
    # Loop through files
    n <- 1
    dtInts <- c()
    dtNames <- c()
    for (rsFile in rsFilelist) {
        # Check dates
        fileStr <- sub(inPath, rsFile, replacement="")
        dtStr <- substr(unlist(strsplit(fileStr,"[.]"))[2], 
                        nchar(unlist(strsplit(fileStr,"[.]"))[2])-(7-1), 
                        nchar(unlist(strsplit(fileStr,"[.]"))[2]))
        #dtPOSIXct <- as.POSIXct(dtStr, format="%Y%j", tz="UTC")
        #begPOSIXct <- if(is.null(begin)) { dtPOSIXct } else { as.POSIXct(begin, format="%Y.%m.%d", tz="UTC") }
        #endPOSIXct <- if(is.null(end)) { dtPOSIXct } else { as.POSIXct(end, format="%Y.%m.%d", tz="UTC") }
        #if ( (dtPOSIXct >= begPOSIXct) & (dtPOSIXct <= endPOSIXct) ) {
        rsRast <- raster::raster(rsFile)
        # Remove MODIS nodata values
        if (!is.null(noData)) {
            if (noDataQual == "min") {
                rsRast[rsRast[] < noData]<-NA
            } else if (noDataQual == "max") {
                rsRast[rsRast[] > noData]<-NA
            } else {
                rsRast[rsRast[] == noData]<-NA
            }
        }
        # Apply MODIS scaling
        if ((valScale != 1) | (valAdd != 0)) {
            rsRast <- rsRast * valScale + valAdd
        }
        # Track dates
        dtInts[n] <- as.integer(difftime(as.POSIXct(dtStr, format="%Y%j", tz="UTC"), 
                                         as.POSIXct("198001", format="%Y%j", tz="UTC", 
                                                    units="days")))
        dtNames[n] <- paste0("DT", format(as.POSIXct(dtStr, format="%Y%j", tz="UTC"), 
                                          "%Y.%m.%d"))
        if (n==1) {
            rsStack <- raster::stack(rsRast)
        } else {
            rsStack <- raster::addLayer(rsStack, rsRast)
        }
        n <- n+1
        #   } # end date check
    } # end for loop
    names(rsStack) <- dtNames
    if (!is.null(outFile)) {
        ConvertStack2NC(rsStack, outFile, varName, varUnit, varLong, varNA)
        return(rsStack)
    }
    return(rsStack)
}


#' Convert a raster stack to a NetCDF file.
#' 
#' \code{ConvertStack2NC} takes a raster stack of RS images and outputs a NetCDF
#' file with a time dimension and specified variable name & metadata.
#' 
#' \code{ConvertStack2NC} converts a raster stack to an output NetCDF file. The
#' raster stack should be already processed through the
#' \code{\link{ConvertRS2Stack}} tool or follow the same layer (date) naming
#' convention.
#' 
#' @param inStack The name of the raster stack to export.
#' @param outFile Name for an output NetCDF file.
#' @param varName Name for the NetCDF export variable.
#' @param varUnit Units for the NetCDF export variable.
#' @param varLong Long name for the NetCDF export variable.
#' @param varNA Value to set for "NA" or "no data". Default is -1.e+36.
#' @param flipY Reverse y indices, e.g., to match geogrid S->N 
#' orientation (DEFAULT=TRUE).
#' @return NULL
#'   
#' @examples
#' ## Export the raster stack of LAI images created through ConvertRS2Stack
#' ## to a NetCDF file. Use the full time series of images.
#' 
#' \dontrun{
#' ConvertStack2NC(lai.b, outFile="BCNED_LAI.nc", varName="LAI", 
#'                 varUnit="(m^2)/(m^2)", varLong="Leaf area index")
#' }
#' @keywords IO
#' @concept MODIS dataMgmt
#' @family MODIS
#' @export
ConvertStack2NC <- function(inStack, outFile=NULL, varName=NULL, varUnit=NULL, 
                            varLong=NULL, varNA=-1.e+36, flipY=TRUE) {
    # Get dates
    dtInts <- c()
    dtNames <- names(inStack)
    for (i in 1:length(dtNames)) {
        dtStr <- sub("DT", dtNames[i], replacement="")
        dtInts[i] <- as.integer(difftime(as.POSIXct(dtStr, format="%Y.%m.%d", tz="UTC"), 
                                         as.POSIXct("198001", format="%Y%j", tz="UTC", 
                                                    units="days")))
        }
    # Output NetCDF file
    raster::writeRaster(inStack, outFile, "CDF", overwrite=TRUE,
            varname=varName, varunit=varUnit, longname=varLong,
            xname="west_east", yname="south_north", zname="Time", 
            zunit="days since 1980-01-01", 
            bylayer=FALSE, NAflag=varNA)
    # Set the time variable
    ncFile <- ncdf4::nc_open(outFile, write=TRUE)
    ncdf4::ncvar_put(ncFile, "Time", dtInts)
    ncdf4::nc_close(ncFile)
    if (flipY) {
      system(paste0('ncpdq -O -a -south_north ', outFile, ' ', outFile))
    }
}


#' Run MODIS-R Whittaker smoothing over pre-processed raster stack.
#' 
#' \code{SmoothStack} takes a raster stack of RS images and outputs a smoothed 
#' raster brick over the same time period.
#' 
#' \code{SmoothStack} converts a raster stack of RS images (as processed through
#' \code{\link{ConvertRS2Stack}}) and calls the MODIS-R
#' \code{\link[MODIS]{whittaker.raster}} smoothing function to generate a
#' smoothed raster brick over the same time series as the input. All function
#' parameters are per the \code{\link[MODIS]{whittaker.raster}} function except 
#' we force the timeInfo to be derived from the input raster stack so the names
#' match and the smoothed brick can be used in other tools. The
#' \code{\link[MODIS]{whittaker.raster}} tool also exports a set of smoothed
#' TIFs, so also specify an output file directory.
#' 
#' @param inStack The name of the raster stack to smooth.
#' @param w FROM \code{\link[MODIS]{whittaker.raster}}: In case of MODIS
#'   composite the 'VI_Quality' raster-Brick, Stack or filenames. Use preStack
#'   functionality to ensure the right input.
#' @param t FROM \code{\link[MODIS]{whittaker.raster}}: In case of MODIS
#'   composite the 'composite_day_of_the_year' raster-Brick, Stack or filenames.
#'   Use preStack functionality to ensure the right input.
#' @param lambda FROM \code{\link[MODIS]{whittaker.raster}}: _Yearly_ lambda
#'   value passed to ?ptw:::wit2. If set as character (i.e. lambda="600"), it is
#'   not adapted to the time serie length but used as a fixed value (see
#'   details). High values = stiff/rigid spline
#' @param nIter FROM \code{\link[MODIS]{whittaker.raster}}: Number of iteration
#'   for the upper envelope fitting.
#' @param outputAs FROM \code{\link[MODIS]{whittaker.raster}}: Character.
#'   Organisation of output files: single each date one RasterLayer, yearly a
#'   RasterBrick for each year, one one RasterBrick for the entire time-serie.
#' @param collapse FROM \code{\link[MODIS]{whittaker.raster}}: logical, if TRUE
#'   the input data is treated as _1_single_Year_ collapsing the data using the
#'   Julian date information without the year.
#' @param outDirPath FROM \code{\link[MODIS]{whittaker.raster}}: Output path
#'   default is the current directory.
#' @param removeOutlier FROM \code{\link[MODIS]{whittaker.raster}}: Logical. See
#'   details
#' @param outlierThreshold FROM \code{\link[MODIS]{whittaker.raster}}: Numerical
#'   in the same unit as vi, used for outliers removal. See details
#' @param mergeDoyFun FROM \code{\link[MODIS]{whittaker.raster}}: Especially
#'   when using argument collapse=TRUE, multiple measurements for one day can be
#'   present, here you can choose how those values are merged to one single
#'   value: "max" use the highest value, "mean" or "weighted.mean" use the mean
#'   if no weighting "w" is available and weighted.mean if it is.
#' @param ... nFROM \code{\link[MODIS]{whittaker.raster}}: Arguments passed to
#'   ?writeRaster (except filename is automatic), NAflag, datatype,
#'   overwrite,...
#' @return raster brick of smoothed images
#'   
#' @examples
#' ## Take the raster stack of LAI images created through ConvertRS2Stack and 
#' ## apply a smoothing filter that also removes outliers, which we specify to 
#' ## be more than 0.5 LAI from the smoothed value.
#' \dontrun{
#' lai.b.sm <- SmoothStack(lai.b, 
#'          outDirPath="/Volumes/d1/adugger/RS/MODIS_ARC/PROCESSED/FRNTRNG_LAI_SMOOTHED", 
#'          outputAs="one", removeOutlier=TRUE, outlierThreshold=0.5, lambda=1000, 
#'          overwrite=TRUE)
#' }
#' @keywords smooth
#' @concept MODIS dataAnalysis
#' @family MODIS
#' @export
SmoothStack <- function(inStack, w=NULL, t=NULL, lambda = 5000, nIter= 3, 
                        outputAs="one", collapse=FALSE, outDirPath = "./",
                        removeOutlier=FALSE, outlierThreshold=NULL, 
                        mergeDoyFun="max", ...) {
    if (!file.exists(outDirPath)) {
      dir.create(outDirPath, showWarnings = FALSE)
    }
    timeInfo <- MODIS::orgTime(inStack, pos1 = 3, pos2 = 13, format = "%Y.%m.%d", 
                               pillow=0)
    resultList <- MODIS::whittaker.raster(vi=inStack, w=w, t=t, 
                                          timeInfo=timeInfo, lambda=lambda, 
                                          nIter=nIter, outputAs=outputAs, 
                                          collapse=collapse, 
                                          outDirPath=outDirPath, 
                                          removeOutlier=removeOutlier, 
                                          outlierThreshold=outlierThreshold, 
                                          mergeDoyFun=mergeDoyFun, ...)
    resultBrick <- resultList[[1]]
    names(resultBrick) <- names(inStack)
    resultBrick
    }


#' Inserts pre-processed images into appropriate forcing NetCDF files by date.
#' 
#' \code{InsertRS} takes a raster stack or brick of RS images and exports
#' individual images to matching (by date) forcing NetCDF files.
#' 
#' \code{InsertRS} takes a raster stack or brick  (as created by
#' \code{\link{ConvertRS2Stack}} or \code{\link{SmoothStack}}) or a NetCDF file
#' (as created by \code{\link{ConvertStack2NC}}) of RS images and exports each
#' layer (time step) to the appropriate time step forcing file. Only looks for
#' the date (not time) and inserts at the 00:00 hour on that date. The input
#' stack, brick, or file should be already processed through the
#' \code{\link{ConvertRS2Stack}}, \code{\link{SmoothStack}}, or
#' \code{\link{ConvertStack2NC}} tools or follow the same layer (date) naming
#' convention.
#' 
#' @param inFile The name of the raster stack/brick or NetCDF file (full
#'   pathname) to export.
#' @param forcPath Path to the forcing data you want to modify. Forcing data
#'   files MUST match the size/resolution of the images in the inFile.
#' @param forcName The suffix for the forcing data files to modify
#'   (DEFAULT="LDASIN_DOMAIN1")
#' @param varName Name for the NetCDF variable to export. The varibale will be
#'   copied as-is, so make sure it matches the variable name needed in the
#'   forcing data.
#' @param varUnit Units for the NetCDF export variable. Only required if the
#'   inFile is a raster stack/brick. If the inFile is a NetCDF file, the units
#'   will carry over.
#' @param varLong Long name for the NetCDF export variable. Only required if the
#'   inFile is a raster stack/brick. If the inFile is a NetCDF file, the
#'   longname will carry over.
#' @param varNA Value to set for "NA" or "no data". Default is -1.e+36.
#' @param overwrite Boolean to allow the tool to overwrite existing variables if
#'   found in the forcing data. (DEFAULT=TRUE)
#' @return NULL
#'   
#' @examples
#' ## Export the raster stack of LAI images created through ConvertRS2Stack to 
#' ## the forcing data.
#' 
#' \dontrun{
#' InsertRS(lai.b, forcPath="FORCING", forcName="LDASIN_DOMAIN3",
#'          varName="LAI", varUnit="(m^2)/(m^2)", varLong="Leaf area index")
#' 
#' ## Export the NetCDF of LAI images created through ConvertStack2NC to the 
#' ## forcing data.
#' 
#' InsertRS("BCNED_LAI.nc", forcPath="FORCING", forcName="LDASIN_DOMAIN3", 
#'          varName="LAI")
#' }
#' @keywords IO
#' @concept MODIS dataMgmt
#' @family MODIS
#' @export
InsertRS <- function(inFile, forcPath, forcName="LDASIN_DOMAIN1",
                        varName=NULL, varUnit=NULL, varLong=NULL, varNA=-1.e+36,
                        overwrite=TRUE) {
    # Raster stack/brick case
    if (inherits(inFile, "RasterStack") | inherits(inFile, "RasterBrick")) {
        dtNames <- names(inFile)
        for (i in dtNames) {
            dtStr <- sub("DT", i, replacement="")
            dtStrForc <- paste0(unlist(strsplit(dtStr,"[.]"))[1], 
                                unlist(strsplit(dtStr,"[.]"))[2], 
                                unlist(strsplit(dtStr,"[.]"))[3], "00")
            ncFile <- ncdf4::nc_open(paste0(forcPath,"/",dtStrForc,".",forcName), 
                                     write=TRUE)
            dimT <- ncdf4::ncdim_def( "Time", "", 1, unlim=TRUE, create_dimvar=T)
            dimY <- ncdf4::ncdim_def( "south_north", "", 1:dim(inFile)[1], create_dimvar=T)
            dimX <- ncdf4::ncdim_def( "west_east", "", 1:dim(inFile)[2], create_dimvar=T)
            # NOTE: ncdf4 reads dimensions in reverse order from ncdump!
            varNew <- ncdf4::ncvar_def(name=varName, units=varUnit, 
                                       dim=list(dimX, dimY, dimT), 
                                       missval=varNA, longname=varLong)
            if ( (varName %in% names(ncFile$var)) ) {
                if (overwrite) {
                    ncdf4::nc_close(ncFile)
                    system(paste0('ncks -O -x -v ', varName, ' ', 
                                  paste0(forcPath,"/",dtStrForc,".",forcName), ' ', 
                                  paste0(forcPath,"/",dtStrForc,".",forcName)))
                    ncFile <- ncdf4::nc_open(paste0(forcPath,"/",dtStrForc,".",forcName), 
                                             write=TRUE)
                } else {
                    stop(paste0('Error: Variable ', varName, ' exists but overwite is set to FALSE. Exiting.'))
                    }
                }
            ncdf4::ncvar_add(ncFile, varNew)
            # Close and re-open, otherwise ncdf4 throws error
            ncdf4::nc_close(ncFile)
            ncFile <- ncdf4::nc_open(paste0(forcPath,"/",dtStrForc,".",forcName), 
                                     write=TRUE)
            ncdf4::ncvar_put(ncFile, varNew, RotateCw(raster::as.matrix(inFile[[i]])))
            ncdf4::nc_close(ncFile)
            }
    # NetCDF file case
    } else if ( inherits(inFile, "character") & file.exists(inFile) ) {
        inNC <- ncdf4::nc_open(inFile)
        nTime <- inNC$var[[varName]]$dim[[3]]$len
        for (i in 1:(nTime)) {
            dtNum <- inNC$var[[varName]]$dim[[3]]$val[[i]]
            dtStr <- as.POSIXlt("198001", format = "%Y%j", tz = "UTC")
            dtStr$mday <- dtStr$mday + dtNum
            dtStrForc <- paste0(format(dtStr, "%Y%m%d"), "00")
            if (overwrite) {
                system(paste0("ncks -A -v ", varName, " -d Time,", i-1, ",", i-1, " ", 
                              inFile , " ", paste0(forcPath,"/",dtStrForc,".",forcName)))
            } else {
                stop('Error: The NetCDF insert option uses ncks append which does NOT check for duplicate variables. Existing variables will automatically be overwritten. To OK this option, please set overwrite=TRUE.')
                }
            }
        ncdf4::nc_open(inFile)
        }
}


#' Calculates summary statistics from a remote sensing time series
#' 
#' \code{CalcStatsRS} takes a raster stack/brick of RS images and generates a
#' time series dataframe of summary statistics across the domain for each time
#' step.
#' 
#' \code{CalcStatsRS} takes a raster stack or brick of remote sensing images
#' over a time period (as created by \code{\link{ConvertRS2Stack}} or
#' \code{\link{SmoothStack}}) and generates a dataframe object that summarizes
#' all cells in the RS image at each time step. The statistics calculated are
#' mean, min, max, and standard deviation. This tool is useful for evaluating
#' how a smoothing function is impacting the time series of images.
#' 
#' @param inStack The name of the raster stack or brick to calculate statistics
#'   on.
#' @return A dataframe of statistics by time period (date).
#'   
#' @examples
#' ## Calculate domain statistics and plot the mean.
#' 
#' \dontrun{
#' stats.lai.b <- CalcStatsRS(lai.b)
#' with(stats.lai.b, plot(POSIXct, mean, typ='l'))
#' }
#' @keywords univar
#' @concept MODIS dataAnalysis
#' @family MODIS
#' @export
CalcStatsRS <- function(inStack) {
    minDf <- as.data.frame(raster::cellStats(inStack, stat="min", na.rm=TRUE))
    maxDf <- as.data.frame(raster::cellStats(inStack, stat="max", na.rm=TRUE))
    meanDf <- as.data.frame(raster::cellStats(inStack, stat="mean", na.rm=TRUE))
    sdDf <- as.data.frame(raster::cellStats(inStack, stat="sd", na.rm=TRUE))
    statDf <- cbind(meanDf, minDf, maxDf, sdDf)
    statDf$POSIXct <- as.POSIXct("1980-01-01", format="%Y-%m-%d", tz="UTC")
    for (i in 1:nrow(statDf)) {
	statDf$POSIXct[i] <- as.POSIXct(sub("DT",row.names(statDf)[i], replacement=""), 
                                  format="%Y.%m.%d", tz="UTC")
        }
    colnames(statDf) <- c("mean","min","max","sd","POSIXct")
    rownames(statDf) <- NULL
    statDf
    }


#' Gap fill a series of MODIS images using interpolation.
#' 
#' \code{GapFillRS} takes a set of pre-processed MODIS TIF files and
#' fills no-data gaps using interpolation.
#' 
#' \code{GapFillRS} scans the specified directory and runs the GDAL python
#' script gdal_fillnodata.py. If a mask file is provided, GapFillRS will
#' convert specified mask cells to specified values AFTER the gap-filling.
#' 
#' @param inPath The path to the directory that holds the already processed
#'   MODIS TIF files.
#' @param matchStr The regular expression for filename matching (e.g.,
#'   "*Lai_1km.tif").
#' @param outDir Pathname to the directory to store the output files.
#' @param begin Date string for the start date to include. The date string
#'   should follow \href{http://r-forge.r-project.org/projects/modis/}{MODIS}
#'   package convention (e.g., "2011.06.01"). (DEFAULT=NULL, all files are
#'   processed).
#' @param end Date string for the end date to include. The date string should
#'   follow \href{http://r-forge.r-project.org/projects/modis/}{MODIS} package
#'   convention (e.g., "2011.06.01"). (DEFAULT=NULL, all files are processed).
#' @param maskFile Mask grid file. File must be in a GDAL-friendly format and 
#'   should match the dimensions and resolution of the inputs.
#' @param maskVals List of value sets to pull from mask file and convert. 
#'   The list should be in the form: list("maskValue"=outputValue), 
#'   e.g., list("1"=0, "16"=0) to convert all cells in the mask grid with mask 
#'   values of 1 and 16 to 0 in the output files.
#' @param nodataVal Value to use to represent no data. Make sure this value 
#'   is NOT in your mask raster. (DEFAULT=-9999) 
#' @param scriptPath OPTIONAL path to where gdal_calc.py script lives in case
#'   it is not in the same default path as gdal executables
#' @param pos1 From MODIS orgTime: Start position of date in the filename (DEFAULT=10).
#' @param pos2 From MODIS orgTime: End position of date in the filename (DEFAULT=16).
#' @param format From MODIS orgTime: How is the date formatted in the file, default 
#' expects: YYYYDDD (\%Y\%j). Read ?as.Date for for more information.
#' @return null
#'   
#' @examples
#' ## Gap fill a set of MODIS LAI images. We will convert all non-vegetation 
#' ## land use types into 0.
#' 
#' \dontrun{ 
#' maskVals <- list("1"=0, "16"=0, "19"=0, "23"=0, 
#'                  "24"=0, "25"=0, "26"=0, "27"=0)
#' GapFillRS("~/wrfHydroTestCases/MODIS_ARC/PROCESSED/FOURMILE_LAI",
#'                          outDir="~/wrfHydroTestCases/MODIS_ARC/PROCESSED/FOURMILE_LAI_GAPFILL",
#'                          matchStr="*Lai_1km.tif", 
#'                          begin="2013.05.01", end="2013.07.31", 
#'                          maskFile="~/wrfHydroTestCases/Fourmile_Creek/DOMAIN/geo_LUINDEX.tif",
#'                          maskVals=maskVals)
#' }
#' @keywords IO
#' @concept MODIS dataMgmt
#' @family MODIS
#' @export
GapFillRS <- function(inPath, matchStr, outDir, 
                            begin=NULL, end=NULL,
                            maskFile=NULL, maskVals=NULL, 
                            nodataVal=-9999, scriptPath=NULL,
                            pos1=10, pos2=16, format="%Y%j") {
  # Setup
  if (!is.null(scriptPath)) scriptPath <- paste0(scriptPath, "/")
  inPath <- path.expand(inPath)
  outDir <- path.expand(outDir)
  if (!is.null(maskFile)) maskFile <- path.expand(maskFile)
  if (!file.exists(outDir)) {
    dir.create(outDir, showWarnings = FALSE)
  }
  
  # Get file list
  if (!is.null(begin) & !is.null(end)) {
    beginDt <- format(as.POSIXct(begin, format="%Y.%m.%d", tz="UTC"), "%Y%j", tz="UTC")
    endDt <- format(as.POSIXct(end, format="%Y.%m.%d", tz="UTC"), "%Y%j", tz="UTC")
    timeInfo <- MODIS::orgTime(MODIS::preStack(path=inPath, pattern=matchStr), 
                               begin=beginDt, end=endDt, pillow=0,
                               pos1=pos1, pos2=pos2, format=format)
  } else {
    timeInfo <- MODIS::orgTime(MODIS::preStack(path=inPath, pattern=matchStr), pillow=0,
                               pos1=pos1, pos2=pos2, format=format)
  }
  rsFilelist <- MODIS::preStack(path=inPath, pattern=matchStr, timeInfo=timeInfo)
  #rsFilelist <- list.files(path=inPath, pattern=glob2rx(matchStr), full.names=TRUE)

  # Mask case
  if (!is.null(maskFile)) {
    # Setup mask calc string
    calcStr1 <- ''
    calcStr2 <- ''
    for (i in 1:length(maskVals)) {
      maskIn <- names(maskVals)[i]
      maskOut <- maskVals[[maskIn]]
      if (calcStr1=='') {
        calcStr1 <- paste0(maskOut, "*(A==", maskIn, ")")
        calcStr2 <- paste0("+(", nodataVal, ")*((A!=", maskIn)
      } else {
        calcStr1 <- paste0(calcStr1, "+", maskOut, "*(A==", maskIn, ")")
        calcStr2 <- paste0(calcStr2, ")&(A!=", maskIn)
      }
    }
    calcStr <- paste0(calcStr1, calcStr2, "))")
    # Create mask
    tmpMask <- tempfile(fileext=".tif")
    cmd_calc <- paste0(scriptPath,"gdal_calc.py ", 
                     "-A ", "'", maskFile, "'",
                     " --outfile=", "'", tmpMask, "'",
                     " --calc='", calcStr, 
                     "' --NoDataValue=", as.character(nodataVal), 
                     " --overwrite")
    system(cmd_calc)
  }
  
  # Loop through files
  for (rsFile in rsFilelist) {
    cmd_fill <- paste0(scriptPath,"gdal_fillnodata.py ", 
                      rsFile, " ", outDir, "/", basename(rsFile))
    system(cmd_fill)
    # Mask case
    if (!is.null(maskFile)) {
      cmd_merge <- paste0(scriptPath,"gdal_merge.py ", 
                       " -o ", "'", paste0(outDir, "/", basename(rsFile)), "'",
                       " -n ", as.character(nodataVal),
                       " '", paste0(outDir, "/", basename(rsFile)), 
                       "' '", tmpMask, "'")
      system(cmd_merge)
    }
  }
  file.remove(tmpMask)
}